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Abstract 

Using the specific model of a bilayer of classical charged particles (bilayer Wigner crystal), 
we compare the predictions for energies and pair distribution functions obtained by 
Monte Carlo simulations using three different methods available to treat the long range 
Coulomb interactions in systems periodic in two directions but bound in the third one. 
The three methods compared are: the Ewald method for quasi-two dimensional systems 
[D.E. Parry, Surf. Sci. 49, 433 (1975); ibid., 54, 195 (1976)], the Hautman-Klein method 
[J. Hautman and M.L. Klein, Mol. Phys. 75, 379 (1992)] and the Lekner summations 
method [J. Lekner, Physica A176, 485 (1991)]. All of the three method studied in this 
paper may be applied to any quasi-two dimensional systems, including those having not 
the specific symmetry of slab systems. For the particular system used in this work, the 
Ewald method for quasi-two dimensional systems is exact and may be implemented with 
efficiency; results obtained with the other two methods are systematically compared to 
results found with the Ewald method. General recommendations to implement with 
accuracy, but not always with efficiency, the Lekner summations technique in Monte 
Carlo algorithms are given. 
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I. Introduction. 



In Chemical physics, Condensed matter physics and Molecular physics, numerical sim- 
ulations are basic theoretical tools to study generic model systems and to predict qual- 
itatively and/or quantitatively the properties of some materials. Molecular simulations 
allow to study the thermodynamical and the transport properties of some more or less 
detailled models of real substances. In Molecular simulations, there are mainly two pro- 
cedures used: the Molecular Dynamics method and Monte Carlo algorithms, both being 
generally considered as complementary [1, 2]. 

The number of particles or the number of degree of freedom of a system generally con- 
sidered in these studies is ranging from few hundreds to few millions. Compared to 
macroscopic materials, this number is small and far from the thermodynamical limit. To 
avoid irrelevant surface or finite-size effects, periodic boundary conditions are generally 
used. 

When the interaction is short ranged, truncations of potentials or minimum image con- 
ventions are used to compute with accuracy the interaction energy or forces between 
molecules or particles. 

On the other hand, when the interaction is long-ranged, as it is the case for Coulombic 
and Dipolar interactions, truncation of the potential may lead to severe bias [3]. For 
bulk- like systems with Coulombic or Dipolar interactions, Ewald methods are employed 
to take into account the long-ranged contribution of the potential to energy or forces 
when periodic boundary conditions are applied to the system [4-6]. 

In the context of molecular simulations, the interesting feature of Ewald summations 
method for bulk-like systems is that the cost in computing time for the calculation of the 
energy of a configuration may scales as N2 , where N is the number of particles. For large 
systems, this computational time is still too expensive and some numerical procedures 
based on the Fast Fourier Transform algorithm like the particle-particle/particle-mesh 
(PPPM) method [7], or procedures based on multipole expansion like the fast multipole 
methods [8] have been implemented. For bulk-like systems, all these methods have been 
used and conveniently intercompared in the last two decades. 

Some very interesting systems may not be considered as bulk-like systems, for instance 
when one or two dimensions of the system are small compared to the other dimensions: 
the quasi-one dimensional or quasi-two dimensional systems. Examples of such systems 
are fluid-fluid or fluid-solid interfaces, monolayers, biological membranes, free standing 
films, cylindrical pores, carbon nanotubes, etc. As for bulk-like systems, to avoid some 
irrelevant finite size effects, other than the finite extension of these systems, partially pe- 
riodic boundary conditions are employed. Partially periodic boundary conditions consist 
of taking periodic boundary conditions for only one or two dimensions, while the remain- 
ing dimensions are considered with their finite extension. For these systems, straightfor- 
ward applications of the Ewald summations lead to a computational effort that scales 
as N 2 and not as iV 3 / 2 , because the analytical formulation of the electrostatic energy of 
this kind of systems is more complicated than for bulk-like systems [9-15]. To compute 
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coulombic or dipolar energies in these quasi-two (or quasi-one) dimensional systems, sev- 
eral others methods have been proposed [9,16-23]. One has to note that some of these 
methods may be applied only to systems having a slab-like geometry, while the methods 
compared in this paper may be applied to any quasi-two dimensional systems. 
The purposes of the present work is to exhibit a comparison between three different 
methods proposed to study quasi-two dimensional systems. The method compared in 
this work are: (a) the Ewald method for quasi-two dimensional systems (Ewald quasi- 
2D) [10,11,13-15], (b) the Hautman-Klein method [16] and (c) the Lekner summations 
technique [19]. All the three methods are applied to a very simple system: a bilayer 
Wigner crystal, described in the section II of this paper. As it is explained further in 
the paper, we have chosen to apply these methods to this particular system because for 
this system the Ewald quasi-2D method may be implemented with the same efficiency as 
the Ewald method for bulk-like systems. To make consistent intercomparisons between 
methods, one needs to have some reliable and solid reference points, the Ewald quasi-2D 
method is the best candidate to produce such reference points. All the results obtained 
by applying the Hautman-Klein method and the Lekner summations technique to this 
system are systematically compared to the results obtained by applying the Ewald quasi- 
2D method. 

This work follows a previous study on the convergence of Lekner summations done by 
the author in ref.[24] and another study done in ref.[25] where we have compared the 
methods of references [16,18,20] on a system having a slab geometry. One of the aims 
of the present work is also to present explicitly a procedure to implement with accuracy 
(but not always with efficiency) the Lekner summations technique in Monte-Carlo sim- 
ulations using the Metropolis algorithms. 

Comparisons of Ewald summation techniques for confined and quasi-bidimensional sys- 
tems have been done previously, such as in the interesting work by A.H. Widmann and 
D.B. Adolf [26]. In these studies, authors generally compare numerical accuracy and 
efficiency of methods on only few configurations of a simple system; for instance in the 
study of ref.[26], Widmann and Adolf compared Hautmann-Klein, Nijboer-de Wette and 
Parry methods on five configurations of a system with 100 particles. In the present work, 
the Ewald quasi-2D, Hautmann-Klein and Lekner methods are implemented in Monte 
Carlo simulations; the influence of the accuracy on thermodynamical averages and on 
correlation functions is outlined. The purpose of the present paper is to concretely illus- 
trate, on a very simple system, how some methods, and especially the Lekner method, 
may introduce biais on physical quantities in a Monte Carlo simulation, since relations 
between thermodynamical averages and accuracy of methods are needed to validate out- 
puts of numerical simulations. 

In section II, we describe the model on which comparisons between the different methods 
are made. In section III, we present the three methods that are studied and compared 
in this paper. Section IV is devoted to exhibit the results obtained with Monte Carlo 
simulations using the Metropolis algorithm. A general discussion and some recommen- 
dations to implement correctly and with accuracy the Lekner summations are given in 
section V. 
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For completeness, an appendix that gives the computation of the interaction between 
two charged surfaces and ions has been added at the end of the paper. 

II. Description of the Model and simulations pa- 
rameters. 

A bilayer of charged particles is the 'reference' system chosen in this work. The three 
methods used to compute the electrostatic energy of this system are presented in the 
next section. 

This model is frequently used to give crude representations of very different systems. For 
instance such bilayer Wigner crystals are used to give representation of strongly coupled 
electronic bilayers of charged particles in two dimensional semiconductor heterostructures 
or in dusty plasmas. Some recent experiments on Laser-cooled 9 Be + ions have provided 
direct observations of structural phase transition in these systems [27] . 
This very simple model is also used to give very crude representation of neutral lipid 
bilayers [28] or as basic theoretical model to facilitate the conceptual understanding of 
counterion-mediated attractions between similarly charged planes [29]. The monolayer 
version of this model is also useful to describe classical electrons trapped on the surface 
of liquid helium [30] . 

Our reference system is a bilayer made of N = 2Nq point ions interacting by a Coulomb 
potential 1/r. The ions are evenly distributed in two layers Li and L2 separated by a 
distance h. A snapshot of an instantaneous configuration of a system of 512 point ions is 
represented on Figure 1. Each layer is a square of side L and partial periodic boundary 
conditions, with a spatial periodicity L, are applied in both directions parallel to the 
layer (directions x and y) while no periodic boundary conditions are taken in the third 
direction (direction z). The purpose of this work is not to study the phase diagram of this 
system (which may be very rich, see for instance ref.[31]), thus in all computations that 
we have done, the shape and area of the layers are constants; nevertheless most of the 
computations done in this paper can be compared to a more extensive thermodynamical 
study of this system done in ref.[32]. The charge of points ions is q and to guaranteed 
charge neutrality, both layers Li and L2 carry a uniform surface charge density a given 
by 

Nq + 2aL 2 = (1) 

In all computations done, we have taken q = 14. The characteristic length a, the ion-disk 
radius, is defined by irpoa 2 = 1 where the ion density in each layer is po = Nq/L 2 . This 
ion-disk picture is a simplified version of the hexagonal Wigner-Seitz cell. 
We have made computations using the three differents methods described in the next 
section in six differents situations resumed in TABLE I. For Runs a, c, d, e, e' and f, the 
coupling constant is T = q 2 /kTa = 196 (a = 1), while for Runs b, it is F = 98 (a = 2). 
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For all computations, a Monte Carlo cycle (MC-cycle) is made of a random trial move 
of each particle in both layers, the trial moves are accepted or rejected according to the 
Metropolis algorithm applied to the NVT ensemble [33, 1, 2]; in a MC-cycle there are 
N trial moves. The amplitude of the trial displacement has been chosen such as the 
acceptance ratios ranged between 30% and 60%. No exchange of ions between layers Li 
and L2 is allowed. From an initial configuration, t eq MC-cycles are used for equilibration 
and after t av MC-cycles are used to accumulate the thermodynamical averages. The 
differences between Runs e and e' are in the values of t eq and t av ; Runs e' were applied 
only to computations using the Ewald quasi-2D and Hautman-Klein methods. 
In Runs a, c, d, e, e' and f, the thermodynamical stable state is a square or hexagonal 
solid; for these computations, the average electrostatic energies per particle are closed 
to the corresponding Mandelung energies. The Mandelung energy, fiuo, for a square 
bidimensional lattice may be found in a computation done by Totsuji [5]. Taking the 
lattice length as L/y/No, the Mandelung energies computed in the work of Totsuji are 
given in TABLE I for all Runs done in the present work. Since the shape and the area 
of the bilayer are kept constant, this would impose preferentially a square lattice for the 
solid-like phase, even if the thermodynamical stable phase has another symmetry. To 
study rigourously the solid phase diagram of this bilayer system, one should allow at least 
the shape of the box to fluctuate. This study was done in ref.[32]. In the present work, 
the shape, a square of side L, and the area L 2 of layers are constant and the average 
energy per particle found in Monte Carlo computations, except for Runs b, is closed to 
uq. We also give in TABLE I, the preferential structures found in the computations done 
in ref.[32] for the bilayer Wigner crystal. 

The order parameters \l/ m used in ref.[32] to differentiate the square and hexagonal 
structures are not computed in the present study. In the spirit of the current work, 
these parameters are irrelevant and since the fixed square shape of layers would induced 
preferentially a square lattice for the solid-like phases, the values for ^4 and ^6 that 
would have been found might be misleading to further studies of the solid phase diagram 
of this bilayer system. 

The energy of the bilayer system can be written as 

U = ^intra ~\~ Winter (2) 

where Ei n t ra is the intralayer contribution and E{ n t er is the interlayer contribution, in- 
cluding the contribution of the surface-surface energy (see appendix). In the next section, 
we give the analytical expressions for Ei ntra and Ei nter for each of the methods studied 
in this work. 

The intralayer gu and interlayer gi% pair distribution functions, corresponding to par- 
ticles in the same layer and to particles belonging to differents layers but positions 
projected onto the same plane respectively, are evaluated. These distribution functions 
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are computed as 
5ii (a) 



L 2 

ieLijeLij^i 



(3) 



L 2 

912 (a) = ^ < S S *( s " s w) > 
ieLi jeL 2 



where particles positions are noted by = Sj + Zje 2 and < . . . > denotes statistical 
averages. 

In the next section, we describe the different methods used to compute the electrostatic 
energy, while in section IV we present the results obtained by applying these methods 
to the bilayer system. 



III. Numerical methods to compute the electro- 
static energy. 

An efficient technique for computing the long ranged part of the intermolecular interac- 
tions for systems with periodic boundary conditions is provided by Ewald sums [1, 2, 4]. 
In this procedure, the interaction energy between molecules is separated into short ranged 
and long ranged contributions by using a damping function f(r; a), where a is a conver- 
gence parameter conveniently chosen. For the coulombic interaction, this separation is 
done by setting 

1 = 1 - f{r;a) + f(r;a) 

The short range contributions to the energy are given by the first term of the right 
handed side (r.h.s) of Eq.(l) and the long ranged contributions are given by the second 
term. To handle the long ranged contributions, Fourier transform of the second term is 
taken; the interaction energy can be set in the form 

E = E ( r s) + E ( £ (5) 

where E^ is the short ranged contribution and E^ the long ranged contribution, which 
for technical convenience, is expressed as 

The damping function f(r;a) may be chosen such as both contributions of Eq.(5) are 
rapidly convergent. Such choices allow to derive efficient algorithm to compute the 
interaction energy. 
The choice 
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2 rocr 

f(r;a) = erf(or) = —= / dt exp(— t 2 ) 
V7T Jo 

is frequently used. Other derivations using an integral representation of the Gamma 
function give the same damping function [13, 15]. We found for the Coulomb interaction 
energy of iV particles in a simulation box of side L with periodic boundary conditions in 
the three directions, 

4 s) = ^EEE^ erfcHrii + Lnl) 



E 



(I) 



i=i j=i n 

N 



\rij + L n \ 



N 

a nt- 2 



1=1 v 1=1 



( _ 2n ^ 



exp(— k 2 /4a 2 



A' 



lE^exp^fc.r^p 

i=i 



(7) 



with 



^x-^x^x ^yLy&y ~\~ TlzLz&z 



where n x , n y and n z are relatives integers. 

The primed sum in E$ indicates that for L n = 0, i = j must be omitted. In Eq.(7), 
Ti = (xi,yi,Zi) are positions of particles, = r« — r 3 - and fc the lattice vectors in the 
tridimensional Fourier space. 

In Eqs.(7), the long ranged contributions E^} q ar e expressed with a summation over the 
particles, N contributions to the sum, and not as a summation over the pair of particles, 
N(N — 1)/2 contributions; this factorization is one of the technical features that allows to 
implement efficiently the Ewald-3D method. By reference to the crystalographic origin 
of the Ewald summations, p{k) = J2j Qj exp(ik.rj) is sometimes called the structure 
factor (of the simulation box) . 

The same procedure, with the same damping function, has been applied to systems in 
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two dimensions [5, 6, 34]. In particular for Coulomb interaction, it gives 

N N I 



with 



i=l j = l V 

N 



eifc(a\sij + L u \) 
\sij + Li/ 1 



E 



(0 

K =0 a A 



— N N 

t(2^0 --j=2^Qi 



i=i 



(8) 



K^O i=l j=l 



-t'l' — v x L x e x + VyL y e.y 



where ^ x and v y are relative integers. 

k are the lattice vectors in the bidimensional Fourier space, Sj = the position of 

the particles in the two dimensional space and = s^ — Sj. The convergence parameter 



a is usually chosen such as the summations that give E r can be restricted to v 
This choice is often taken as aL ~ 5 — 8. 

For systems with a quasi-two dimensional geometry, as the bilayer Wigner crystal de- 
scribed in section II, the Ewald method for three dimensional systems, given by Eqs.(7), 
may be applied by taking a highly asymmetric box [35, 20, 21(b)]. For systems with a 
slab geometry, this procedure is certainly the most efficient [2 1(b), 24]. A careful exami- 
nation of the summations may show that large errors, arising by a naive implementation 
of the tridimensional Ewald method can be corrected analytically by adding some cor- 
rections terms [20,21 (b)]. Nevertheless, these procedures may not easily be applied to 
systems that do not have the specific symmetry of slab-like systems, as it is the case 
for fluid-fluid, fluid-solid interfaces, monolayers, etc. For these reasons, several methods 
have been developed to take into account the finite extension of systems without any 
periodicity in the direction of the finite extension. 

The method for real two dimensional systems, given by Eqs.(8), may of course not be 
applied to quasi-two dimensional systems, except as a limiting case when the length of 
the finite extension in the third direction tends to zero (h — > 0). 



0. 



A. The Ewald method for quasi-two dimensional systems. 

Following analytical computations done in refs. [10-15], one may derived an Ewald method 
for three dimensional systems with periodicity in only two directions. 
One of the most simple derivation of the Ewald quasi-2D method is the analytical deriva- 
tion done by Parry [10] in 1975; this derivation of the Ewald quasi-2D method is also 
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very useful to link this method to the Ewald method for tridimensional systems. In the 
derivation done by Parry, the Ewald-3D method is taken as a starting point. By letting 
the periodicity of the periodic boundary condition in the direction where the system has 
a finite extension tending to infinity, an analytical formulation of the Ewald quasi-2D 
method is found. Other analytical derivations, that do not use the Ewald-3D method 
as a starting point, may also be found in refs. [11,13,15,35]. In particular, an elegant 
method may be found in a recent work by Grzybowski, Gwozdz, and Brodka [15]; their 
derivation is also closely related to the rapidly convergent representations for Green's 
functions for Laplace's equations presented in a work by Linton [37]. 
We define our notations as following, 













= L u + n z L z e z 


k 


— & ~~\~ k%€iz 



where k z = 2tt/L z . 

One finds by taking L z — > oo 

E 2D = E r + E K=Q + E^q = Eff (9) 

with the contributions given by [10] 

1 erfc(a\rij + L u \) 
E r = 2 E E E W \r-- + L,A ( M 

z i= i j= i v \n 3 + ^v\ 



vr N N / exp(— a 2 z 2 -) \ a N 

E k=0 = ~\ E E 9iQj( ^ 11 + l%|erf(a|% |)) - -y= E 9? ( 10 - b ) 

i=i j=i » » i=i 

N N /■ x 

= ^EEE wn*, « ; Zij f-^pd (10 . C ) 

K^0*=ii=i 

where 

F(k, a; Zij) = exp(KZjj)erfc(^ + aZij) + exp(-KZjj)erfc(^- - azij) (11) 
and the contribution k = arising from fe z / 0. 

When applied to the bilayer Wigner crystal described in section II, one may use the fac- 
torization as one particle summation in the long ranged contributions given by Eq.(lO.c). 
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Thus, one may implement with the same efficiency the Ewald quasi-2D method for this 
system as it is done in the Ewald-3D method. This situation is very particular to the 
system considered in the present work. For other quasi-2D systems, because the com- 
plicated form of the function F(k, a; Zij), one may not in general find a factorization in 
one particle summation as in the Ewald-3D method. For the bilayer Wigner crystal, Zij 
may take only two values: if both particles are in the same layer, then = or else 
Zij — h. 

The intralayer contributions to the long ranged part of the electrostatic energy (z^ = 0), 
for k / are given by 

E kSo 4E erfC(K/2a) ( I E *exp(iK. S< ) I 2 + I E ft exp(iK.-<) I ") (12) 

and for k = 

= K £ + ( S "'ft - 7? { S + £ ?? ) (13> 

Since the electroneutrality of the system is achieved with the surface charge density 
of layers Li and L2, contributions with (X^ g L 2 ft) 2 7^ have to be included in the 
electrostatic energy (see appendix for more details). 

The interlayer contributions to the long ranged part of the electrostatic energy for k / 
are given by 

E k% = T E »[( E ftexp(iK. Si ))( E (14) 

k^O ieLi jeL 2 

where F(k, a; h) is given by Eq.(8) and $l[z] is the real part of the complex number z. 
For k = 0, the interlayer contribution given by Eq.(lO.b) is 

= -£( + E «)( E *)■ (") 

ieLi jeL 2 

The electrostatic energy of a configuration of the bilayer Wigner crystal with the Ewald 
quasi-2D method is given by 

171 t?2D 1 771 / Tjiintra , ^intra , T7iintra\ , / winter , winter , ciinter , 771 \ 
h = h + h w = {h r + h K=Q + h K¥ jQ ) + {h r + h K=Q + + h w ) 

with Eqs. (lO.a), (12-15) and (A.14). 

This expression of the electrostatic energy is an exact application of the Ewald method 
to the bilayer system; all the computations with Monte Carlo algorithm done with this 
method are considered in section IV as reference points. 

Generally, one may not implement with efficiency the Ewald quasi-2D method, thus there 
are very few uses of this method. Of course, the procedure described in this subsection 
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has been used in the study of the solid-phase diagram of the bilayer crystal [32] and also 
in other numerical studies of electronic bilayer systems [38]. Most applications of the 
Ewald quasi-2D method in computer simulations use precalculated tables of potential 
energy, forces and second derivatives on a three-dimensional grid and calculations are 
performed by interpolation of the tables [39, 35, 20]. 

In ref.[40], Arnold and Holm have shown that the MMM2D method has a computation 
effort that scales as iV 5 / 3 (log N) 2 . It is interesting to note that for the bilayer Wigner 
system, it is possible to implement the Ewald quasi-2D method with the same efficiency 
as the Ewald 3D method, the computation effort scaling as iV 3 / 2 . Thus, for this system, 
the Ewald quasi-2D method is slightly more efficient than the MMM2D method. 



B. The Hautman-Klein method. 

The computation of particle-particle interaction in the Hautman-Klein method starts 
with the identity [16] 

I I M z 2n M z 2n 

n=0 n=0 

with 

(-l) n (2n)! 



a. 



2 2 "(n!) 2 



By using damping functions h n (s;a) on terms l/s 2n+l , one separates the electrostatic 
energy of the system in a short ranged part E s and long ranged part E\ , 



E HK = E S + E l = Eff 



where 



1 N N ' / 1 M hi \\ 

e s = \ E E w E - E -nzfr^m^) (17) 
ei = \ E E w E «n4 ra ( E K{ ^r ) (is) 

z i=l j=l n=0 {V} ijM 

As in the Ewald method, the long ranged contributions are evaluated by using a Fourier 
transform. In the original derivation by J. Hautman and M.L. Klein the damping func- 
tions h n (s;a) are chosen such as 

ho(s; a) = erf(as) 
s 2 ^ 1 a n (2n)\ K s > 



and 
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The derivation done by J. Hautman and M.L. Klein with the damping functions given in 
Eqs.(19) allows to recover easily the analytical derivation of the Ewald method for exact 
two dimensional systems given in Eqs.(8) by letting Zij — > 0. This is obtained obviously 
from equations (17-19). Using the damping functions of Eqs.(19), we have 



1 N (I 



erf(as 



1 



i=1 j¥=i 



M 

S( 2 ")' 



v 2n 



v 2n ( 



erf(as 



v 



■)) 



(20) 



and 



N N M 

A i=lj=l 



Jin 



^ 2n)ru £ K 2 «- 1 erfc(K/2a)^ 



in. Si 



(21) 



v i=i 



According Eqs.(16) for n = 1 to 3 we have 



hi(s; a) = erf(as) — 
h 2 (s; a) 
hs(s; a) = erf(as) — 



2as 

. . 2as 2 2 
erf(as) —e ° 



'(l + 2aV) 



2as 



(1 + 2 a 2 s 2 _ 4 4 4 8 6 6 

v 3 9 9 ; 



(22) 



225 + 225 jj 

The attractive feature of the method is that, by writing the zfj 1 explicitly as polynomials 
in Zi and Zj, the long ranged contribution given by Eq.(21) can be cast into a sum of 
terms each of which involving product of two functions having the general form 



N 



F p( K ) = ^2li z i exp(iK.Sj) 



(23) 



A careful examination of Eqs. (20-21) in one hand and Eqs.(10) in the other hand, shows 
that Hautman-Klein method may not be considered exactly as Taylor expansion of the 
Ewald quasi-2D method. In particular, the complicated contribution E k= q in Eq.(lO.b) 
in the Ewald quasi-2D method is roughly approximated in the Hautman-Klein method 
by a constant. As a consequence, the contributions to intralayer and interlayer energies 
do not have exactly the same constant terms. 

For the bilayer Wigner crystal, the long ranged intralayer contributions E 1 ^^ HK are 
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exactly given by Eq.(12) and with Eq.(21), one may expressed ^^q hk as Eq.(13). 
Long ranged interlayer contributions are given by 

E ^6hk = ~ A E *[( E ftexp(iK. Si ))( E «,-exp(-«. 8i ))] (24) 

/c^U ieLi jel^2 

with the same notations as in Eq.(14) and 

M -J 

a; fc) = ]T — - («fc) 2n erfc(K/2a) (25) 

n =0 \ in )- 

The contribution to interlayer energy for k = 0, according to Eqs.(21) and (13) has to 
be written as 

41 e oV = -^f(E^)(E^)- (26) 

instead of the Taylor expansion of Eq.(15). 

The electrostatic energy of a configuration of the bilayer Wigner crystal with the Hautman- 
Klein method is given by 

TTi jpHK , i7 / rnintra , zriintra , j-iintra \ 

h = h + h w - ^ r HK + & k=Q hk + b K ^ HK ) 

i / pintcr _j_ pinter , pintcr , 771 \ 
+ K^r,HK + ^ K =0,HK + h K^Q,HK + h W) 

with Eqs.(20), (12), (13), (24), (26) and (A.14). 

Because an expansion of 1/r for small z is taken as a starting point in Hautman-Klein 
method, this method is often considered as inaccurate for rather thick systems. Despite 
its obvious inaccuracy, for slab-like systems, as the system studied in ref.[25], the results 
found with the Hautman-Klein method are in very good agreement with results obtained 
with a purely electrostatic model and results obtained with the methods of refs. [18,20]. 
The Hautman-Klein method may also be implemented with almost the same efficiency 
as Ewald methods for bulk-like systems, but without taking any periodic images in the 
direction of the finite extension of the system. Consequently, the Hautman-Klein method 
is specifically well adapted to systems as fluid-fluid or fluid-solid interfaces, that do not 
have a slab-like geometry. The Hautman-Klein method may also be implemented to 
take into account electrostatic images at the interface between two media with different 
dielectric constants. This method has been used recently in a study of a mixed lipid/non- 
ionic surfactant membrane [41]. 



C. The Lekner method. 

The third procedure used in the present work is a method proposed by J. Lekner [19]. 
This method may also be more or less easily related to the original Ewald method. 
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When applied to quasi-one dimensional systems, Lekner summations method is an eigen- 
function expansion of the Green function of Laplace equation with periodicity in only 
one direction (see for instance Eq.(3.10) in ref.[37]). Following the derivations done in 
refs. [13, 15,37], an Ewald method for quasi-one dimensional systems may be derived [42], 
the Lekner summations for the quasi-one dimensional systems may be recovered by let- 
ting the convergence parameter a, defined in Eq.(4), tends to infinity [37, 43]. 
For quasi- two dimensional systems, the relation between Lekner summations and Ewald 
summations is sligthly more complicated. With the derivation that permits in quasi- 
one dimensional systems to relate the Lekner summations to an eigenfunction expansion 
of Green functions, one is obtaining for quasi-two dimensional systems the so-called 
Nijboer-de Wette representation [12, 13, 44], the same method as the one used recently 
in the method called MMM2D (see ref.[21(a)] and [45]); Lekner summations are related 
to these methods by an asymmetric use of the Poisson-Jacobi formula. As for quasi-one 
dimensional geometries, when the convergence parameter a tends to infinity, one may 
relate easily the Nijboer-de Wette representation (or the 'far formula' of the MMM2D 
method [21]) to the Ewald method [37]. 

Following the original derivation done by J. Lekner [19] and computations done by N. 
Gr0nbech- Jensen and co-workers [46-49] and by S. Marshall [50], the electrostatic energy 
of a quasi-two dimensional system can be computed as 

1 N N N 

i=ij=ij^i i=i 

The interaction energy Ey between pairs is approximated in a numerical computation 
by one of the two formulas 



Vij(n c ,n K ) 



or 

Uij(n c ,n K ) 




The integers n c and nx are the truncation parameters that would govern the accuracy of 
the finite summations [24, 37]. The main technical problem that one would encounter in 
trying to implement straightforwardly this procedure is in the convergence rate of one of 
the two summations for some special configurations of the pair of particles i, j relatively 
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to the geometry of the simulation box. For these particular configurations, that are very 
frequent when few billions of configurations are sampled, one of the two summations 
given by Eqs.(28) and (29) is very slowly convergent and a lot of contributions must be 
included to compute the correct value of the energy. The origin of this slow convergence 
rate is in the behavior of the modified Bessel function Kq as its argument tends to zero 
(see ref.[24]). To implement correctly the Lekner summations in a Monte Carlo algorithm 
or in a Molecular dynamics simulation, one must use a procedure to overcome this poor 
convergence rate or else some bias would plague the simulation. To achieve this aim, 
several methods have been proposed. A procedure frequently used consist in applying 
a cyclic symmetry that is found in the original derivation of the method [19]. In the 
following, we restrict ourself to this procedure and the implementation of the Lekner 
summations according to this procedure is called the Lekner-cyclic method. At the end 
of this subsection we will review briefly other procedures proposed to overcome the poor 
convergence rate of one of the summations Uij(n c ,riK) or V^(n c , n^)- 
The cyclic symmetry consists in recognizing that as n c — ► oo and nx — ► oo, we have 

lim Vij(n c ,n K ) = lim Uij(n c ,n K ) = Eij (30) 

For the particular configurations where one of the two summations is slowly convergent, 
to compute Eij one may use Vij(n c ,nK) or Uij(n c ,riK) (refs. [19,47,50]). Since for given 
n c and nx, Vij(n c ,riK) an d Uij(n c ,riK) may have very different numerical values, the 
application of the cyclic symmetry may introduce complicated bias [24]. In particular, 
if n c and uk are badly chosen, the value for found in computing energy of particle % 
in the electrostatic potential of the particle j may be different of the value found for Eji 
when the energy is computed by considering the particle j. An inequality as E^ ^ Eji 
happens here because the configuration of the vector Sij, relatively to the simulation 
box, is different from the configuration of the vector Sji, especially when the minimum 
image convention is used. 

The study done in ref.[24] is helpful to define a criterion on s^j to choose to compute 
Vij(n c ,nK) rather than Uij(n c ,riK) or the converse. 

For the bilayer Wigner crystal, the interaction energy E^ between pair of particles is 
computed as 

Vij(n c ,n K ) if \xij\ > \yij\ 

(31) 

Uij(n c ,n K ) if \xij\ < \yij\ 



E^ 



The criterion used here is quite simple because shapes of the layers are squares of side 
L. 

General formulas as Eqs.(28) and (29) for Lekner summations, when the base of unit 
cell is not a square (L x / L y ) or do not have a rectangular shape, have been proposed 
in refs. [48,51-53]; for these geometries a careful examination of the criterion given by 
Eq.(31) should be done, especially when the base of the unit cell is not rectangular. 
To compute, with the cyclic symmetry, the energy of a configuration made of N particles, 
the number of different contributions to evaluate is of the order n c (2n^-|-l)A 2 . As it was 
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shown in our previous work [24], when the criterion given by Eq.(31) is used, for some 
configurations of the vector Sy the value computed by using the most rapidly convergent 
of the two formulas (28) and (29) becomes independent of the direction of the vector sij 
for rather high values of n c (e.g. n c ~ 40 for a relative accuracy of 1CT 5 — 1CT 6 ) [24, 37]. 
On the other hand, for most of the configurations of the vectors Sij, the convergence of 
summations (28) and (29) are more rapids, thus a value as large as n c ~ 40 is not needed 
in all computations done with the Lekner summations. Then, to improve the efficiency 
of the method, n c is not chosen independently of the argument of the Bessel function. 
Taking into account the asymptotic behavior of the Bessel function Kq(x) 



K o( x ) - \/J ex P(" 



-x) as x — ► oo 



truncation to nx = 3 or 4 is enough to achieve a good accuracy (.Kg (19.) ~ 1.6 10 
The value of n c used is chosen such as 



1 < n c (sij) < n™ ax (32) 

where n™ ax is a prefixed value and n c (Sy ) the value of n c chosen in such a way that the 
output of the summation, Eq.(28) or Eq.(29), would be independent of the vector Sij. 
In the computations presented in the next section, n c (sij) is chosen as the first integer 
such as 

2nn c (s ll )[{ I f + k) 2 + {f) 2 Y 2 >l9- (33) 

If the integer found by using Eq.(33) is greater than n™ a:r , then one takes n c {sij) = n^ iax . 
For each value of the term between the brakets in Eq.(33), a value of n c (sij) is taken 
from a precomputed table that takes into account the prescription imposed by n™ ax . 
This evaluation of n c (sij) is done for all pairs of particles. For the bilayer Wigner crystal 
this procedure is used for both the intralayer and interlayer contributions. 
This procedure has been applied to all computations done with the Lekner-cyclic method 
presented in section IV, in this section n c denotes the value of n™ aa: used in Eqs.(32) 
and (33). 

For the bilayer Wigner crystal, the separation of the energy into intralayer and interlayer 
contributions is done as in the Ewald quasi-2D and Hautman-Klein methods as 

E™ = E^ + Ej£ e[ = Efi>) (34) 

with 



« a = ^E E ^E E ^+E^ elt +E^ elt (^) 



where 



i 2 2 00 00 

Uf di = - r lim ( Eij - ^) = | [4 E E K (^) + 7 - ln(4vr)] (36) 



L 

m=l k=l 
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and 

Winter = E E (37) 
For the bilayer Wigner crystal, we have 

vi^(n c ,n K ) =4^f>s(^^) £ i^m)^ + *|] 

m=l k=-n K ^gg^ 

-^ln[l-cos(2^)]-Ml ln2 
for intralayer contributions, and 

^ ter K,n*) = 4^f>s(^) g KoM(f + fe) 2 + (^)f] 

m=l '" k=—riK 

-ML In [cosh (2tt£) - cos (2t3)] - ^ In 2 

(39) 

for interlayer contributions and similar relations for C/|j ltra and ?7™ ter (interchanging 
with j/y). 

The electrostatic energy of a configuration of the bilayer Wigner crystal with the Lekner- 
cyclic method is given by 

E = E LC + E w = ^tra + Winter + E w) (40) 
with equations (35) to (39) and (A. 14). 

The Lekner summations method has been used in few studies using molecular simulation 
algorithms. 

The cyclic symmetry has been used recently in a Molecular Dynamics study of thin 
water-acetonitrile films [54]: in this work an appendix that describes the method used to 
implement the Lekner-cyclic method may be found. In recent Monte Carlo simulations 
of adsorption of colloidal particles on a charged surface [55], the Lekner-cyclic method 
is also used [56], but technical details on the implementation are not given in ref.[55]. 
To implement correctly the Lekner summations method, other procedures have been 
proposed. In the present work, none of these procedures, reviewed briefly below, have 
been neither implemented, nor compared to the Lekner-cyclic or Ewald quasi-2D meth- 
ods. 

In ref.[57], J.F. Harper proposed to use an integral representation of summations when 
arguments of Bessel functions are small. This method might be well adapted for quasi- 
one dimensional systems [42] where, instead of having two different analytical relations, 
Eqs.(28) and (29), one has only one summation. In these situations one cannot use the 
cyclic symmetry. 

Another procedure is to derive an alternative expression which converges fast as the ar- 
gument of the Bessel functions tends to zero. Using the Hurwitz zeta function, R. Sperb 
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had proposed such alternative expressions [58]. The Sperb method has been used in a 
Molecular Dynamics study of strongly charged polyelectrolyte brushes [57]; in this work 
one would find an appendix that describes the procedures used by authors. In ref.[43], A. 
Grzybowski and A. Brodka also used a procedure close to the Sperb method for quasi-one 
dimensional systems. The Lekner-Sperb method was also used by A.G. Moreira and R.R. 
Netz in their studies of counterions distributions near charged plates [60, 61]. In ref.[61], 
Moreira and Netz give a detailled and very clear description of their implementation of 
Lekner summations. In their implementation they use Sperb transformation with the 
following procedure: if p = + k) 2 + z' 2 > 1/3, they use relations similar to Eqs.(28) 
and (29) of the present work with truncation parameters (n c ,nx) — (11,2), if p < 1/3 
Sperb formula is used. As it will be shown in section IV, when Lekner-cyclic method 
is used, truncation parameters as (n c , nx) — (10,3) do not give accurate results. Nev- 
ertheless, one should not conclude that the truncation parameters used by Moreira and 
Netz are too small to give accurate results; in the core of their method they use Sperb 
formula which main purpose is to reduce these truncation parameters. Implementation 
done in ref.[61] and in the present work are not exactly identicals. 

Another procedure has been developed by J.N. Newman based on two-dimensional 
Chebyshev expansions and economized polynomials [62]. 

There are few others studies that used Lekner summations on very different, interesting 
and complicated systems [55,63-67], but details on implementation are not given there. 
Although the cyclic symmetry is certainly the procedure the mostly used for quasi-two 
dimensional systems, details on the implementation are needed. Firstly, there are at 
least three different methods to implement the Lekner summations method. Secondly, 
as it will be shown in section IV, the results of a computation may be depending on 
parameters and criterions used to correct the slow convergence rate of one of the two 
summations given by Eqs.(28) and (29). Moreover, the criterions used to implement 
the Lekner cyclic method in Eqs.(31), (32) and (33) are depending on the geometrical 
parameters of the simulation box or on parameters of the studied system. 

IV. Results. 

In this section, we present the results obtained with the three methods of the previous 
section. 

The three computer codes used in this work are very similar. The main differences be- 
tween the three programs are in the subroutines that compute the electrostatic energy, all 
others parts of the programs are the same. In particular, the files in which configurations 
are written may be used in any of the three programs, the random numbers generator 
used is the same and it is possible, with our codes, to start computations for each of 
the three methods from exactly the same initial condition and with exactly the same 
sequence of random numbers for the trial moves of the particles and for the acceptation 
tests. This point was necessary not only to compare the compatibility of the methods, 
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but also to test the efficiency and the accuracy of the methods. 

This section is made of two subsections. In subsection IV. A, using the same initial config- 
uration and the same random numbers sequence, we compare the sampling of the phase 
space of the bilayer Wigner crystal done with the three methods. In subsection IV. B, the 
results of computations for the Runs defined in TABLE I are given. All the computations 
were performed on a Silicon Graphics Origin 2000 (R10000/195 processors). 

A. Accuracy of the sampling of the phase space. 

In the Monte Carlo Metropolis algorithm [1, 2, 33], the phase space of the system studied 
is sampled by an importance-weighted random walk. From an initial configuration of 
the system, its phase space is sampled by a construction of configurations in such a way 
that configurations are generated with a probability proportional to the statistical weight 
of the configurations. A convenient way is to generate a trial configuration from an old 
configuration of the system and to accept or reject the trial configuration according to an 
acceptance probability defined with the help of the detailed balance condition. According 
to this procedure, if the trial configuration is accepted, it becomes the new configuration 
of the system, or else, if the trial configuration is rejected, the old configuration is 
taken as the new configuration. The procedure is implemented by building another trial 
configuration from the new configuration. 

When the trial configuration is generated only from an old configuration and when the 
acceptance probability may only be defined with the knowledge of the old and trial 
configurations, the chain of configurations generated by the Metropolis algorithm has 
all the properties of a Markov chain. All ensembles of the statistical mechanics can be 
sampled with the Metropolis algorithm [1]. 

To apply the Monte Carlo algorithm to the canonical ensemble, the simplest procedure 
is to generate trial configurations from old configurations by a random displacement of 
one particle of the system. The displacement of the particle (the new configuration of 
the system) is then accepted according to a probability given by 

acc(o — > n) = min(l, exp(— [3(E n — E )) (41) 

where E n is the total energy of the trial (new) configuration and E Q the total energy 
of the old configuration. For these procedures, the accuracy of the computation of the 
energy is closely related to the accuracy of the sampling of the phase space of the system. 
For the bilayer Wigner crystal, the energy difference used in the acceptance probability of 
Eq.(41) is computed by using one of the three methods of section III and in the following 
we set 

0AEW = -P(E n - E ) (42) 

Since for the Lekner-cyclic method, configurations where the point ions are on a regular 
two dimensional lattices are very penalizing for the convergence of summations, these 
kinds of configurations were not taken as initial configurations. To generate initial con- 
figurations on which an examination of the influence of the truncation parameters n c 
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and rix could be done, we have chosen to build initial configurations as follows. N$ 
point ions are arranged on a square lattice in each layers; this configuration is taken 
as the initial configuration for the code using the Ewald quasi-2D method. From this 
initial configuration, 100 MC-cycles were performed with the Ewald quasi-2D method. 
After these 100 MC-cycles (200iVo trial moves), the point ions are no longer arranged in 
a perfect regular lattice, thus this configuration may be taken as an initial configuration 
for all the three methods. This procedure presents also the advantage of allowing to 
start computations from a well defined configuration that we may reproduce easily, if 
the same sequence of random numbers is taken. 

If codes are written such as the random numbers are generated following the same se- 
quence and used for the same purpose in each of the three codes, then, if the three 
methods give exactly the same energy, the trajectory in the phase space obtained by the 
Markov chain of configurations would be exactly the same in the three computations. 
The trajectory followed by the system with the Ewald quasi-2D method is considered 
as the exact trajectory; the trajectories obtained with the Hautman-Klein and Lekner- 
cyclic methods are compared to this trajectory. In TABLE II, we give some elements 
to appreciate the accuracy of the sampling of the phase space by the Hautman-Klein 
and Lekner-cyclic methods relatively to the sampling obtained with the Ewald quasi-2D 
method for initial configuration that corresponds to Runs a (cf. TABLE I). The energy 
difference (3AEp P for the first trial move of the first MC-cycle from the same initial 
condition with the same sequence of random numbers is given for different truncation 
parameters for the Hautman-Klein and Lekner-cyclic methods. This first move is partic- 
ularly interesting, whatever it is accepted or rejected, the move is identical in the three 
methods and the position of all the other particles are also the same in the three meth- 
ods. The values of (3AEpp v ^ show that the truncation parameter n c in the Lekner-cyclic 
method governs the accuracy of the computations, while nx seems to have less incidence 
on accuracy. In agreement with our previous work [24], to reach a relative accuracy of 
10 -4 — 10~ 5 , one needs to choose n c ~ 40, while uk = 2 or 3 seems to be sufficient. A 
careful examination of other energy differences shows that, for few particular moves, one 
needs to take = 3 to avoid a lost of accuracy; nevertheless these events are rare and 
in almost all cases, uk = 2 might be sufficient. On the other hand, n c < 25 always leads 
to important differences in almost all energy differences. 

Because energy differences in the three methods are not exactly the same, the trajectory 
in the phase space followed by the system are not exactly the same. In TABLE II, we 
give also the number of trial moves accepted after the first MC-cycle (t=l) and after 
the first hundred MC-cycles (t=100) from the same initial condition. In a deterministic 
point of view, the trajectory in the phase space, using the Metropolis algorithm, has 
to be independing of the method used to compute the energy. Of course, in almost all 
molecular simulations, such deterministic point of view is unreachable, at least because 
of round-off errors done by the machines on floating points [68] . In our study, as soon as 
a random trial move in the computation using the Ewald quasi-2D method is accepted 
(or rejected), while it is rejected (or accepted) in computation using Hautman-Klein of 
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Lekner-cyclic methods, the trajectories are no longer the same. Since in all situations ex- 
aminated, one may not exactly follow the same trajectory with the Hautman-Klein and 
Lekner-cyclic methods as with the Ewald quasi-2D method, we should find a criterion 
to estimate if these methods are enough accurate to give reliable statistical averages. 
To answer this question, one has to come back to the principles of the Monte Carlo sam- 
pling. To obtain reliable thermodynamical averages by making a Monte Carlo sampling, 
one has to generate configurations with a probability that would be proportional, with 
a rather good accuracy, to the statistical weight of the configurations. Thus, it is not 
necessary that the trajectories followed by the system will be exactly the same, even if 
the initial configurations and the sequence of random numbers are the same. The impor- 
tant point is to generate configurations with a probability proportional to the statistical 
weight of the configurations; if this is not achieved, the thermodynamical averages are 
biased in a more or less complicated way. To estimate the bias from the 'exact' results 
given by the Ewald quasi-2D method, we have chosen to follow the instantaneous en- 
ergy per particle. In Figure 2 (a-d), we have plotted for the first 500 MC-cycles the 
instataneous energy per particle obtained with the Ewald quasi-2D on one hand, and 
with others methods on another hand, for the same initial configuration and with the 
same sequence of random numbers. The intantaneous energy of the initial configuration 
is given at t = 0. 

In Figure 2. a, the energies using Hautman-Klein method (M = 3) and Ewald quasi-2d 
are shown. An examination of the first 50 MC-cycles shows that both trajectories are 
exactly the same for these first moves. After, it happens that a trial move is accepted 
in one of the two methods, while it is rejected in the other. Trajectories separate at 
this point. Nevertheless, fluctuations of this quantity are still in good agreement in the 
two methods. This is an indication, but not a demonstration, that both methods would 
give same statistical averages and thus it seems that there are no bias, relatively to the 
Ewald quasi-2D method, in using Hautman-Klein method for the geometrical parameters 
corresponding to Runs a. 

The situation in Figure 2.b, where the same quantities are plotted for Ewald quasi-2D 
method and Lekner-cyclic method with (n c ,nx) = (10,3), is different. Trajectories sep- 
arate in the first trial moves of the first MC-cycle (see also TABLE II). For the first 50 
MC-cycles, fluctuations of the instantaneous energy seem to be similar, but after 100 
MC-cycles the trajectories clearly separate, much more than two magnitudes of the fluc- 
tuation. In this case, as it will be shown in the next subsection, the statistical averages 
are different. 

Since for n c = 10 and uk = 3 as truncation parameters in the Lekner-cyclic method 
the accuracy is not correct, we made the same numerical analysis for higher value of 
these parameters. This is shown on Figure 2.c and 2.d. In Figure 2.c the truncation 
parameters are n c = 25 and uk = 3; the trajectories do not separate for the few first trial 
moves but only after the few first MC-cycles, despite the fact that the relative accuracy 
on the energy difference is only 10~ 2 compared to the energy difference found for the 
Ewald quasi-2D method. For the 500 MC-cycles shown in this Figure, with n c = 25 and 
riK = 3, trajectories stay in good agreement. If one increases the value of n c , as it is done 
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in Figure 2.d where n c = 50, one would find a better agreement between Ewald quasi-2D 
and Lekner-cyclic methods. For n c = 50 and uk = 4, the separation occurs only after 60 
MC-cycles and after the trajectories stay very close. For both cases, (n c , uk) = (25,3) 
and (n c ,riK) = (50,4) one may hope to obtain the same statistical averages than in the 
Ewald quasi-2D method. Of course a dramatic separation, as for (n c , uk) = (10, 3), may 
still happen for both (n c , uk) = (25, 3) and (n c , nx) = (50, 3) after more MC-cycles, this 
had not happened for the first 10000 MC-cycles examined. 



B. Energy and pair distribution functions. 

To estimate statistical averages, we have mainly considered two differents cases. A case 
where the separation between the two layers is small {h = 1.0) for which the coupling 
between the layers is important: this corresponds to Runs a and b. We have also 
considered another case where the separation between layers is large enough (h = 4.0) 
so that both layers are almost independent [32], Runs c to f. 

Since for h = 1.0, the coupling between the layers is important, this has an incidence on 
the interlayer pair distribution functions g\2- This separation between the layers allows 
us to estimate the agreement between methods for small z. On the other hand, because 
for h = 4.0 the coupling between the two layers is very small, one would find for almost 
a ll s i 9i2( s ) — 1- One may estimate on these computations with h = 4.0 the influence of 
the periodic boundary conditions by varying N. 

In TABLE III, the results for Runs a and b are given, U/N is the average energy per 
particle where U is computed in the three methods as 

U =<E> (43) 

Einter/N and Ei n t r a/N are respectively the interlayer and intralayer average energies 
per particle computed according section III; and o\j is an estimation of the statistical 
fluctuation of the total energy computed as 

al =< E 2 > - < E > 2 (44) 

On Figures 3 to 6, we give the intralayer gu(s) and interlayer 512(5) pair distribution 
functions obtained with the Ewald quasi-2D and other methods for Runs a and b. 
The value of o\j given in TABLE III can be used to give an estimation of the accuracy 
needed on the average energy. The fluctuation ajj is the average of the fluctuations 
observed on Figures 2(a-d). The results given on TABLE III show that the Ewald 
quasi-2D and Hautman-Klein methods in one hand and the Lekner-cyclic method with 
(n c , uk) = (25, 3) on the other hand are in a very good agreement for an estimation 
of the energy for Runs a and b. Using the results of TABLE III, one may check that 
the difference between the average energy obtained using Ewald quasi-2D method and 
Hautman-Klein method or Lekner-cyclic method with (n c , n^) = (25, 3) is lesser than 
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2 <jjj. For Runs a, one has | UEwaid ~ Uhk I— 0.2 ay for Hautman-Klein method 
and | UEwaid — Ui e kner |— 0.4 for Lekner-cyclic method with (n c ,%) = (25,3), 
while for Runs b, one has | UEwaid — Uhk |— 0.4 cr^ for Hautman-Klein method and 
I UEwaid — Uiekner |— 1-7 o\j for Lekner-cyclic method with (n c ,riK) = (25,3). The lost 
of accuracy in the Lekner-cyclic method with (n c , tik) = (25, 3) for Runs b may be well 
understood by observing that as L increases the argument of the Bessel functions in 
Eqs.(38) and (39) decreases, thus one has to increase n c to achieve a better accuracy. As 
shown on Figures 3 and 5 the pair distribution functions are also in good agreement, but 
with a small lost of accuracy for Runs b in both Hautman-Klein method and Lekner- 
cyclic method with (n c , nx) = (25, 3) compared to the Ewald quasi-2D method. With 
n c raised to 40, one obtains a very good agreement between the Ewald quasi-2D and 
Lekner-cyclic method for Runs a and b, as it is shown on Figures 6. 
With the same criterions, it appears that the results for the Lekner cyclic method with 
(n c , tlk) = (10, 3) are not in a good agreement with the Ewald quasi-2D method neither 
for Runs a, nor for Runs b. For Runs a, one has | UEwaid — U Lekner |— 4.2 ay (see also 
Fig.2(b)) and for Runs b | UEwaid ~ & 'lekner | — 5.6 ay. Then, as it was mentioned in 
the previous subsection, the truncation parameters (n c ,n/^) = (10,3) are badly chosen 
to give reliable statistical averages in this implementation of the Lekner-cyclic method. 
The bias induced by a bad choice of the convergence parameters in Lekner summations 
is exhibited on Figs.4(a-b). The pair distribution functions, computed with the Lekner- 
cyclic method (n c , nx) = (10, 3) for Runs a, differ dramatically from the pair distribution 
functions computed with the Ewald quasi-2D method. As it is shown on Figs. 5 (a-b), 
when the value of the truncation parameter n c is raised to 25, pair distribution functions 
gn and g\2 computed with the Lekner-cyclic method are exactly the same as those 
computed with the Ewald quasi-2D method. For Runs b, the observations done on 
Figs.4(c-d) and Figs.5(c-d) are not as clear as it was for Runs a; nevertheless one may 
not conclude that for the thermodynamical point corresponding to Runs b, truncation 
parameters (n c ,nft-) = (10,3) are enough to avoid bias, since the average energy is very 
badly estimated. This point could be misleading when testing the implementation of 
the Lekner method in a computer code where no reference points, as those provided by 
the Ewald quasi-2D method, can be found. As mentioned in TABLE I, for Runs b the 
preferential structure of the layers is a disordered fluid-like structure, thus there is less 
configurations of pairs of particles with Xij ~ or y^j ~ compared to a crystal-like 
ordered structure. As a consequence, summations given by Eqs.(38) and (39) converge 
rapidly for most of pairs of particles and the bias induced by the poor convergence rate of 
summations for few pairs of particles is hidden. This situation becomes very dangerous 
when a fluid-solid equilibrium is studied with these parameters. We must outlined also 
that even if the average energy is badly evaluated with truncation parameter n c = 10, 
the value found with this parameter is not very different from the correct value given by 
the Ewald quasi-2D method, because of the oscillatory behavior of summations given by 
Eqs.(38) and (39) (see reference [24]). This point could also be misleading when tests on 
implementation of the Lekner method are done. 

For h = 4.0 (Runs c to f), the statistical averages are given on TABLE IV with same 
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notations. All these computations are done for the same thermodynamical point but 
with different values of the number of particles ranging from N = 128 to 968. According 
to the results of TABLE IV, the Ewald quasi-2D and Hautman-Klein methods are in a 
good agreement despite the large separation of layers (h = 4.0) and except for Runs c 
with N=128. For Runs d, e' and f, one has | Jj T Ewald — Uhk |< 0.7 ajj. On the other 
hand, for Runs c, where N = 128 we have | U Ewald — Uhk |— 3.5 ay and the value 
found for Ei nte r is not correct with the Hautman-Klein method. This is illustrated on 
Figs. 7 (a-b) where it is shown that the pair distribution functions computed with Ewald 
quasi-2D and Hautman-Klein methods are not in agreement. By using the Hautman- 
Klein method with the geometrical parameters of Runs c, one finds a coupling between 
both layers, while with the Ewald quasi-2D method no coupling appears. This lost of 
accuracy in the Hautman-Klein method for Runs c is introduced by the expansion of 1/r 
for small z, one has h/L ~ 0.3 for N = 128, while for N > 338, we have h/L < 0.2, but 
with the same density of point ions. Figures 7 and TABLE IV show that to increase the 
accuracy of the Hautman-Klein method one may increase N. 

The finite-size effect is also outlined on Figures 7 for the Ewald quasi-2D method. The 
pair distribution function gu is sensitive to the finite size of the simulation box, this is 
illustrated in Fig. 7(a) up to Fig. 7(g) where N is ranging from 128 to 968. 
Comparison between Ewald quasi-2D and the Lekner cyclic methods are done in TABLE 
IV and in Figures 8 and 9, for Runs c to f. As for Runs a, the Lekner-cyclic method with 
(n c , uk) = (25, 3) is in very good agreement with the Ewald quasi-2D method. For Runs 
c to f, one has | U E waid-U Lekner |< 1.5 ay for N > 338 and | U Ewa i d -U Lekner |~ 1.8 ay 
for N = 128. For (n c ,n K ) = (10,3), TABLE IV shows that as N is increased, the 
Lekner-cyclic method lost its accuracy. For Runs c with (n c ,nx) = (10,3), one has 
I U Ewald - U Lekner |~ 1.5 ay and for Runs f one finds | U Ewa ld ~ U Lekner |~ 12 ay. This 
point is well illustrated on the representation of g\\ given on Figures 8. As N is increased 
from 128 to 968, the structure functions computed with the Lekner-cyclic method become 
poorly evaluated. With (n c , tik) = (10,3), for N = 128, gu and g\2 are in very good 
agreement with the functions computed with the Ewald quasi-2D method for N = 128 
(Figs. 8 (a-b)), but as N is increased, this accuracy is lost. This lost of accuracy is similar 
to the one observed in Runs b for the Lekner-cyclic method with (n c , nx) = (25, 3) since 
as N is increased, L is increased and thus the argument of the Bessel functions decreased. 
This point may also introduce bias when an implementation of the Lekner method is 
done. Before making extensive computations, one frequently checks the accuracy of 
the computation on systems having very few particles. For systems with long ranged 
interactions, a good estimation of some kind of Mandelung constants is sometimes taken 
as a test of the method [19, 69]. Thus, with this kind of tests, one may conclude using 
systems with very few particles that a small value of n c gives accurate results, but this 
accuracy would be lost as the number of particle is increased. Tests on systems with a 
small numbers of particles (i.e. N < 150) are not sufficient to check the accuracy of the 
implementation of the Lekner-cyclic method. As it is shown on Figure 9, if one raises n c 
to 25, the Lekner-cyclic and the Ewald quasi-2D methods are in very good agreement 
for Runs c to f, both for the average energy and pair distribution functions. 
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On TABLES III and IV, we give also a crude estimation of the average CPU-time 
in seconds per MC-cycle. These estimations are depending on machines and codes. 
Improvement of the efficiency may certainly be achieved; nevertheless these CPU-times 
may be intercompared to give an estimation of the relative efficiency of the methods. 
The average CPU-time per MC-cycle found for the Ewald quasi-2D method in TABLE 
IV may not be connsidered as a representative efficiency of the method, since for the 
particular system considered in this work the Ewald quasi-2D is as efficient as an Ewald- 
3D method (see subsection III. A). TABLE IV shows that the Lekner-cyclic method 
may become particularly expensive as N is increased (see also ref.[61]). With Eq.(33), 
one may estimate the value of n™ ax to achieve enough accuracy with the Lekner-cyclic 
method as N is increased, one roughly has 

n rnax _ jyi (45) 

2-Kh 

Thus, the computing time with the implementation proposed in subsection III.C scales 
like 

St ~ N a with 2 < a < I (46) 



V. Discussion 



The study done in this paper gives a comparison between Lekner, Ewald quasi-2D and 
Hautmann-Klein methods when implemented in a Monte Carlo simulation. In particular 
the influence of the methods on thermodynamical and structural quantities is outlined 
on a very simple model on which the Ewald quasi-2D method may be efficiently imple- 
mented. Thus, the Ewald quasi-2D method and the bilayer Wigner system may serve as 
reference tools to gain on physical quantities of general interest, concretely (i.e. when 
implemented in Monte Carlo Metropolis or in Molecular Dynamics simulations) the accu- 
racy and efficiency of all other methods proposed to handle quasi-bidimensional systems. 
One of the main results of this paper is to demonstrate that the Ewald quasi-2D method 
and the Lekner summation technique are in complete agreement. In refs. [51,67], authors 
have compared the Lekner method with the Hautman-Klein method; the outputs of the 
computations presented in this paper agree with the results of these authors. 
According the computations done here, we may not conclude that the Lekner summations 
are an efficient alternative to the Ewald quasi-2D method, at least for the implementa- 
tion that uses the cyclic symmetry (Lekner-cyclic method). 

The procedure to implement the Lekner-cyclic method, described in the subsection III.C 
of the present paper, is successful to use correctly the Lekner summations in a Monte- 
Carlo Metropolis algorithm. This procedure was based on a previous work where a study 
of the convergence of the Lekner summations is done [24]. The truncation parameter n c 
defined in ref.[24] and in Eqs.(28) and (29) governs the accuracy of the sampling of the 
phase space of a system when Lekner summations are used in Monte-Carlo Metropolis 
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algorithms. Some severe bias may be induced when the method is straigforwardly imple- 
mented. These bias are generated by the diverging behavior of the Bessel functions as its 
argument tends to zero. We may not recommend to use a precise value of the truncation 
parameter n c to always implement the Lekner method with accuracy, because the value 
of n c that must be used is depending on the geometrical parameters of systems studied. 
Thus, for each system a convergence study is to make before performing intensive com- 
putations. 

As explained in subsection IV. B, some tests may be misleading to check the accuracy 
of the implementation. One of the most difficult behavior to handle is the fact that the 
needed value of n c to achieve good accuracy is depending on the number of particles, 
because the spatial periodicity appears explicitly in the argument of the Bessel function 
(see also Figures 8) . The finite size effects in Lekner summations are much more difficult 
to handle, than in Ewald methods; especially because of the asymmetry of Eqs.(28) and 
(29). 

In Ewald methods, it is required that energies are not depending on the convergence 
parameter a. The same requirement is to achieve with the truncation parameter n c in 
Lekner methods. 

The proposed procedure in this present work to implement the Lekner method in molec- 
ular simulations is accurate, but unfortunatly this procedure is not very efficient. Others 
procedures have been proposed to correctly implement the Lekner summations [55, 56, 
60]. These procedures have not been considered in the present work. Some of these 
procedures may improved the efficiency of the method. 

The bilayer Wigner crystal used as a reference model in this work is very convenient, 
since with the Ewald quasi-2D method exact results may be obtained. Thus, this model 
and the results of TABLE III and TABLE IV may be used as references to test the 
accuracy of other methods for quasi-two dimensional systems. 

Lekner summations have also been used to compute Mandelung constants [70]. In these 
computations, the Lekner summations are useful and may be considered as an alterna- 
tive to Ewald methods, mainly because the spatial periodicities that are used in these 
computations are much smaller than the spatial periodicities used in molecular simula- 
tions. 
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Appendix : Electroneutrality and computation of the 
interaction between two charged surfaces and ions. 



For the system described in Section II of the present work and more generally for a 
system made of N charged particles of charge q and two parallel surfaces (Li U L2) 
separated by a distance h with a surface charge density <r(r), electroneutrality is given 
by 

Nq+ f f dra(r) = (A.l) 
J JLiuL 2 

Using conventional notations for the periodic boundary conditions, the electrostatic en- 
ergy of this system is given by 

' air) 



a 2 N N ( ' 1 \ N r 
z i=i 7 =i n \ r i3^ n \ i=1 JL 1 uL 2 n 



r + n 



+- / drdr' > , , 
2 JLxuLa n r r 



<7(r)<r(r0_ ( A - 2 ) 
+ n I 



= -Spp + -EpS 1 + Ess 



The three contributions to the electrostatic energy are: (a) E pp the particle-particle 
interaction energy evaluated in this work by using one of the three methods (i.e. Ewald 
quasi-2D, Hautman-Klein or Lekner-cyclic methods) where self energies are included 
in Epp, (b) E p s is the particle-surface energy and (c) Ess the surface-surface energy, 
including self energy of both surfaces. 

In the following, we consider only uniform surface charge density and using Eq.(A.l), 
one has 

' = - & < A - 3 > 

With these notations, the particle-surface energy is given by 

,2 N 



_ Nq< 
E P s - --fit 



E/ ^Ei r (A.4) 



Following the derivation done by Grzybowski and co-workers in ref.[15] and using the 
periodic boundary conditions, we have 

E pS (e) = -^yv^E (/.(I Zi-\\) + h{\ Zi + \ I)) (A.5) 
i=i 



where 

dt - r»- z2 
3 

is 



roo Af „ \ p -Z 2 e 

f e (\ Z\)= -e~ z 1 = 2v^F -=+ I Z I erf(v^ | Z |) - 20F | Z | (A.6) 

Je £3 L Ve7T J 
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The prescription e has been introduced to outline the diverging behavior of the integral 
in the second member of Eq.(A.6). As e — > 0, we have 

f £ (\ Z |) = 4 - | Z | +2y^ | Z | 2 +o(e§) (A.7) 

As extensively commented in the literature, because of the electroneutrality of the sys- 
tem, the diverging behavior as e — > 0, is cancelled by adding all contributions to the 
electrostatic energy (including self energy contributions). For systems such as all parti- 
cles are confined between charged surfaces (i.e. —h/2 < z% < h/2), as the system studied 
in this work and the system studied in ref.[25], we have 

*i ~ \ I) + U\ z i + ^\) = ^=-V^h + 2^(| Zl - k -? + \ Zl+ h - | 2 ) + o(ei) (A.8) 

Following the same steps, the surface-surface contribution is 

Ess(e) = + -j^r^Uh) (A.9) 

When contributions Ess{e) an d E p s(e) are added, the diverging behavior is not yet 
cancelled by using electroneutrality since one has to include the contribution of the 
particle-particle interaction with the same prescription. 
With the notations of subsection III. A, one has 



E PP (e) = E r + E^ + q -j^Y: E / + (A-10) 

Computing the integral in the right handed side of Eq.(A.lO), one finds 
E pp (e) = E r + E k _^q + E k= q 

2 „2„ N N , -z? e v (A.ll) 

+T^E E (^+l%|erf(v^|%|))+< / ( e ) 



12 S 3 ^w^ 

The particle-particle contribution is then separated into the converging contribution 
given in Subsection III. A and a diverging contribution as 

E pp {e) = E^ + E^(e) (A.12) 

with 

4*> w -^^ + ^(f; £ 4)^ + o (e f) (A.i3) 

Thus, adding the contributions given by Eqs.(A.5), (A.9), (A.ll) and (A.13), and letting 
e — > 0, Eq.(A.2) becomes 

£ = £7^) + | = + 2^I 2 /i = + ^ (A14) 
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The contributions of the charged surfaces to the energy when the particles are confined 
between surfaces is simply a term that is not depending on the position of the particles. 
During the sampling of the phase space of the system described in section II, especially 
when the energy difference for a trial move of a particle is computed, the term added in 
Eq.(A.14), E\y is constant as long as h, L and ./V are constants. This constant is included 
in the results given in TABLE III and IV, and Figures 2(a-d). Ey/ is considered as a 
contribution to interlayer energies. 

On the other hand, is evaluated using the Ewald quasi-2D, Hautman-Klein or 

Lekner-cyclic methods. 
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List of Tables 



TABLE I: Geometrical and simulation parameters of the different Monte Carlo 
samplings. L is the spatial periodicity, h the distance between the layers, N the 
number of point ions in the simulation, a the uniform surface charge density used to 
achieve electroneutrality, t eq the number of Monte Carlo cycles done in each Runs 
for equilibration, t av the number of Monte Carlo cycles to accumulate thermody- 
namical average, (3uq the Mandelung energy computed for a square bidimensional 
lattice. The preferential structure is the structure found in the Monte carlo study 
of the solid-phase diagram of the classical bilayer crystal done in ref.[31]. 

TABLE II: Estimation of the accuracy of the sampling of the phase space by 
the three different methods. fiAE^ is the energy difference computed with the 
three methods of the first trial move of the first MC-cycle from the same initial 
conditions and with the same random numbers. ACC(t = 1) is the number of 
accepted trial moves after the first MC-cycle and ACC(t = 100) the number of 
accepted trial moves after the first hundred MC-cycles (see also Figures 2). 

TABLE III: Averages energies computed with the Monte Carlo algorithm for 
Runs a and b (h — 1.0). (3U /N is the average total energy per particle, (3Ei nte r/N 
and /3Ei ntra /N are respectively the average interlayer and intralayer energy per 
particle, f3ajj/N an estimation of the statistical fluctuation a\j of the total energy 
and 5t is a crude estimation of the average CPU-time in seconds per MC-cycle 
(P = 1/kT, T temperature). 

TABLE IV: Same as TABLE III, but for Runs c, d, e, e' and f (h = 4.0). 
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List of Figures 



Figure 1: Representation of the bilayer Wigner cristal. Periodic boundary 
conditions with a spatial periodicity L are not represented. This representation is 
a snapshot of an instantaneous configuration of a system of 512 point ions with 
the geometrical parameters of the Runs e and e' given in TABLE I. 

Figure 2: Instantaneous energy per particle for the first 500 Monte-Carlo 
cycles computed with different methods but with exactly the same initial configu- 
ration and with the same sequence of random numbers. The parameters of these 
computations are: h = 1.0, N = 512, L = 28.36 (Run a) (see also TABLE II). (a) 
Comparison of the sampling done with the Ewald quasi-2D and Hautman-Klein 
method, (b) Comparison between the Ewald quasi-2D and Lekner-cyclic (n c = 10, 
riK = 3) methods. 

Figure 2: (c) Comparison between the Ewald quasi-2D and Lekner-cyclic 
(n c = 25, riK = 3) methods, (d) Comparison between the Ewald quasi-2D and 
Lekner-cyclic (n c = 50, n K = 3) methods. 

Figure 3: Representation of the intralayer gn(s) and interlayer gn(s) pair 
distributions functions obtained with the Ewald quasi-2D and Hautman Klein 
methods for Runs a and b (h — 1.0). (a) gu(s) from Runs a; (b) gi2{s) from Runs 
a; (c) <7n(s) from Runs b; (d) gn(s) from Runs b. 

Figure 4: Same as the Figure 3, but for the Ewald quasi-2D and Lekner-cyclic 
(n c = 10, n K = 3) methods for Runs a and b (h — 1.0). (a) gn(s) from Runs a; 
(b) 012(5) from Runs a; (c) gu(s) from Runs b; (d) gu(s) from Runs b. 

Figure 5: Same as the Figure 3, but for the Ewald quasi-2D and Lekner-cyclic 
(n c = 25, riK = 3) methods for Runs a and b (h — 1.0). (a) gn(s) from Runs a; 
(b) 512(5) from Runs a; (c) gu(s) from Runs b; (d) g± 2 (s) from Runs b. 

Figure 6: Same as the Figure 3, but for the Ewald quasi-2D and Lekner-cyclic 
(n c = 40, riK = 2) methods for Runs a and b (h — 1.0). (a) gn(s) from Runs a; 
(b) ^12(5) from Runs a; (c) <?n(s) from Runs b; (d) gn(s) from Runs b. 
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Figure 7: Same as Figure 3, but for the Ewald quasi-2D and Hautman Klein 
methods for Runs c to f (h — 4.0). (a) gn(s) from Runs c; (b) gi 2 (s) from Runs c 
(N = 128); (c) g n (s) from Runs d ; (d) g 12 (s) from Runs d (N = 338); (e) g n (s) 
from Runs e and e'; (f) gi 2 (s) from Runs e and e' (N = 512); (g) gu(s) from Runs 
f ; (h) g 12 (s) from Runs f (N = 968). 

Figure 8: Same as Figure 3, but for the Ewald quasi-2D and Lekner-cyclic 
(n c = 10, fix = 3) methods for Runs c to f {h = 4.0). (a) gu(s) from Runs c; (b) 
gi 2 (s) from Runs c (N = 128); (c) gu(s) from Runs d ; (d) gi 2 (s) from Runs d 
(N = 338); (e) gn(s) from Runs e and e' ; (f) gn(s) from Runs e and e' (N = 512); 
(g) (?n(s) from Runs f ; (h) gi 2 (s) from Runs f (N — 968). 

Figure 9: Same as Figure 3, but for the Ewald quasi-2D and Lekner-cyclic 
(n c = 25, tik = 3) methods for Runs c to f (h = 4.0). (a) gu(s) from Runs c; (b) 
gi 2 (s) from Runs c (N = 128); (c) gn(s) from Runs d ; (d) gi 2 (s) from Runs d 
(N = 338); (e) gn(s) from Runs e and e' ; (f) g\ 2 (s) from Runs e and e' (N = 512); 
(g) gn(s) from Runs f ; (h) gi 2 (s) from Runs f (N — 968). 
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TABLE II . MAZARS, Lekner and Ewald sums. 
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Name L h N a t eq t moy (3uq Preferential 

structure 

a 28.36 1.0 512 -4.456 50000 100000 -215.64794 staggered rhombic 

b 56.72 1.0 512 -1.114 50000 100000 -106.68091 fluid-like, disordered 

c 14.18 4.0 128 -4.456 400000 1000000 -215.64794 staggered hexagonal 

d 23.04 4.0 338 -4.456 100000 200000 -215.64794 staggered hexagonal 

e 28.36 4.0 512 -4.456 50000 100000 -215.64794 staggered hexagonal 

e' 28.36 4.0 512 -4.456 200000 1000000 -215.64794 staggered hexagonal 

f 38.99 4.0 968 -4.456 50000 100000 -215.64794 staggered hexagonal 



TABLE I , MAZARS, Lekner and Ewald sums. 



Method Runs L (3U/N (3E inter /N (3E intra /N pau/N St(s) 



Ewald quasi-2d 


a 


28 
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1.0 
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Lekner-cyclic 
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Lekner-cyclic 
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28 
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(n c = 40, n K = 2) 
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56 
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TABLE III , MAZARS, Lekner and Ewald sums. 
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c 


128 


-215.57 
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Lekner-cyclic 


c 


128 


-215.42 


-0.00(2) 


-215.41(9) 


0.100 


0.3 


(n c = 10, n K = 3) 


d 


338 


-215.39 


-0.00(4) 


-215.38(2) 


0.061 


1.7 




e 


512 


-215.86 


-0.00(4) 


-215.85(6) 


0.048 


3.7 




f 


968 


-216.27 


-0.01(2) 


-216.25(6) 


0.045 


13.5 


Lekner-cyclic 


c 


128 


-215.39 


-0.00(2) 


-215.38(9) 


0.101 


0.3 


(n c = 25, n K = 3) 


d 


338 


-215.61 


-0.00(3) 


-215.60(4) 


0.057 


1.8 




e 


512 


-215.67 


-0.00(4) 


-215.66(9) 


0.046 


4.0 




f 


968 


-215.77 


-0.00(5) 


-215.76(8) 


0.034 


14.6 



TABLE IV . MAZARS, Lekner and Ewald sums. 




Figure 1, MAZARS, Lekner and Ewald sums. 
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Figure 2(a),(b), MAZARS, Lekner and Ewald sums. 
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Figure 2(c),(d), MAZARS, Lekner and Ewald sums. 
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Figure 3 (a-d), MAZARS, Lekner and Ewald sums. 
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Figure 4 (a-d), MAZARS, Lekner and Ewald sums. 
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Figure 5 (a-d), MAZARS, Lekner and Ewald sums. 
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Figure 6 (a-d), MAZARS, Lekner and Ewald sums. 
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